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A new approach to the modening of wetting fronts in porous media on the Darcy scale is developed, based on 
considering the types (modes) of motion the menisci go through on the pore scale. This approach is illustrated 
using a simple model case of imbibition of a viscous incompressible liquid into an isotropic porous matrix with 
two modes of motion for the menisci, the wetting mode and the threshold mode. The latter makes it necessary 
to introduce an essentially new technique of conjugate problems that allows one to link threshold phenomena on 
the pore scale with the motion on the Darcy scale. The developed approach (a) makes room for incorporating 
the actual physics of wetting on the pore scale, (b) brings in the physics associated with pore-scale thresholds, 
which determine when sections of the wetting front will be brought to a halt (pinned), and, importantly, (c) 
provides a regular framework for constructing models of increasing complexity. 



1. Introduction 



Th e dynamics of 



1988; 



Olbricht 



1996; 



wetting fronts in 



Alava et al. 



Dorous media remains the subject of intensive research (jAdler &: Brenner 
20041 ). Its main motivation comes, first of all, from a host of important applica- 
tions, notably in oil recovery, hydrogeology and more recently also in carbon dioxide sequestration, microfluidics 
and fuel cells. This topic also poses some fundamental questions about the modelling of evolutionary processes 
in systems with complex topology. In a practically relevant case where the scales of the pore-level and the 
global flow are well separated, one can use an intermediate scale to introduce averaged macroscopic quantities 
and apply the modelling approach of continuum mechanics. In this case, for an incompressible viscous liquid 
invading an isotropic porous medium the averaged flow velocity u and pressure p, both functions of the position 
vector r and time t, satisfy the equation of motion in the form of Darcy's law 
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where k and fi are the permeabihty of the porous matrix and the fluid's viscosity, respectively, and fl is part 

of the porous medium occupied by the fluid. Darcy's law together with the continuity equation V • u = form 

a closed system adequately describing, on the macroscopic level ('Darcy scale'), the bulk distributions of p and 

u. Combining these two equations, one has that the pressure in the flow domain fl is harmonic 



V^p^o, (ren), 



(1.2) 



so that, to determine p, one has to specify two boundary conditions on the part of the boundary of fl whose 
location is unknown (i.e. the wetting front; hereafter dfli) and one boundary condition on the part whose 
position is known (hereafter dfl2, so that dfli U dfl2 = Oft is the boundary of ft). Hence on the wetting front, 
in addition to the kinematic condition 

df 



dt 



u-V/ = 0, (redn^), 



(1.3) 



which specifies the evolution of dfli in terms of its a priori unknown location /(r, t) = 0, we need to formulate 
an appropriate dynamic boundary condition for p. If the dynamics of the displaced fluid also needs to be 
considered, as, for example, in the case of one viscous liquid displacing another, one will still need a dynamic 
boundary condition relating the two fluids' pressures at the interface. 

The main issue in determining the dynamic boundary condition is to what extent the actual physics of wetting 
on the pore scale is acc ounted for, and how i t is represented, on the Darcy-scale level. In particular, as has been 



known for a long time (|Huh fc Scriven 



1971 : 



Dussan V. fc Davis 



1974 : 



Dussan V, 



19791) . the classical model of 



fluid mechanics does not allow viscous fluids to spread over a solid surface with a contac t angle less than 180°, 



whereas numerous experiments show that they do (e.g. see Ch. 3 of iShikhmurzaevi 



20071 ). The way one chooses 



to overcome this ('moving contact-line') problem for the menisci that collectively form the wetting front in a 
porous medium is one of the factors determining how realistic the resulting model will be. At present, this aspect 
of the wetting front dynamics modelling received almost no attention, even in the approach which considers the 
porous medium as a n etwork of capillaries and hence is potentially capable of capturing exac tly the details of 



the pore-scale physics (jLenormand et al. 



1988; 



Aker fc Malov 



2000; 



Joekar-Niasar et al 



20101) 



An alternative to the topologically transparent sharp-interface approach outlined above is to consider the 
wetting front as a transition zone where the volumetric concentration of the displaced and displacing fluid change, 
with no distinction bet ween continuou s and discrete p hases, and treat this, essentially multiphase, system as 



a multicomponent one f| Richards 



1931 : 



Leverett 



19411) , using a thermodynamic closure to relate the pressure 



difference between the two fluid phases with saturations and, possibly, other variables (JHassanizadeh fc Gray 



1993 : 



Mitkov et al. 



1998 : 



Deinert et al. 



20081) . A difficulty in this approach is that, as experiments aimed at 
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Figure 1. Schematic illustration of the meniscus motion in the wetting mode (a) and one of scenarios associated with 

the threshold mode (b). 

determining the closing thermodynamic relationships indicate the need to bring in more and more parameters of 
state, it is not clear whether or not thermodynamics is an adequate tool to describe this, essentially mechanical, 
system. 

2. The model 

In this paper, we introduce a new approach to formulating models describing the propagation of liquid- 
fluid interfaces across porous media based on considering the types (modes) of motion which the menisci that 
collectively form the free surface go through as they advance across the porous matrix. This approach offers a 
regular way of building models of increasing complexity accounting for the topological and geometric features of 
the porous medium. It is worth noting that Darcy's equation (jl.ip in the bulk is itself essentially a consequence 
of the flow profile being approximately parabolic on the pore scale, i.e. it also can be seen as based on a particular 
type of flow. In this sense, the approach we are developing here considers the boundary conditions conceptually 
in the same way as the bulk equations. 

We will illustrate the new approach using the simple case of a viscous incompressible liquid displacing an 
inviscid dynamically passive gas from an isotropic homogeneous porous medium. On the pore scale, each menis- 
cus intersects the pore boundary at a 'contact line', forming a certain 'contact angle' with the solid. We can 
schematically represent the motion of the meniscus on the pore scale as having two principal modes: (i) the 
wetting mode, where the contact line moves across the pore boundary with negligible variation in the meniscus 
shape, and (ii) the threshold mode, where the contact line becomes pinned whereas the meniscus deforms until 
the contact angle it forms with the solid reaches a critical value at which the contact line can move again (Fig.[T]). 
These two modes control the motion as each of them is capable of bringing individual menisci and hence the 
wetting front as a whole to a halt. In the wetting mode this can happen when the contact angle becomes equal 
to the equilibrium one, so that the meniscus no longer needs to move to reach an equilibrium state, whereas in 
the threshold mode, where the contact line is pinned, the flow stops when the pressure building on the meniscus 
(later referred to as p\dQ,i) is insufficient to break through the threshold. 
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Macroscopically, for the pressure on the wetting front measured with respect to the (presumed constant) 

pressure in the displaced gas one has 

Pldfii == Aipi + A2P2, (2.1) 

where pi, p2 are the averaged pressures and Ai, A2 are the spatio-temporaUy averaged fractions of the unit 
area of the free surface corresponding to the two types of motion {Ai + A2 = 1). Importantly, Ai and A2 do 
not have to remain constant as the wetting front propagates, and they are yet to be specified. 

We will begin by considering the wetting mode. For simplicity, we will represent the pore where the wetting 
motion takes place as having a circular cross-section. Then, for low capillary and pore-scale Bond numbers, the 
meniscus will have the shape of a spherical cap and 



Pi — ^2a cos 9d/ a, 



(2.2) 



where a is the liquid-gas surface tension, a is the effective radius of the pore and 6d is the dynamic contact 
angle that the meniscus forms with the solid wall. As shown by experiments on dynamic wetting, in a general 
case 9d depe nds on, and hence should be regarded as a func tional of, the flow field in the vicinity of the moving 



contact line ([Blake et al. 



1999; 



Clarke fc Stattersfield 



20061) . For flow in a porous medium, the contact angle's 



dependence on the fiow field reduces to its dependence only on the contact-line speed ui, which is the leading 
factor determining the flow field in the vicinity of the contact line. Then, 9^ can be described just as a function 
of the contact-line speed: 

Od - F{ui/Uci), (2.3) 

where Ud is an appropriate scale for the velocity. In principle, this dependence, where, for the wetting mode. 
Ml coincides with the sp eed of the meniscus as a whole, could be determined empirically. The theory of fiows 



with forming interfaces ([Shikhmurzaev 



20071 ) specifies the inverse of F{ui/U) as 






(1 + (1 - pf Jcosg,)(cos^, - cosOdY ^ '^' 
4(cos6is-hB) (cos 6*^-^5) 



(2.4) 



where B = {I ~ p\^) ^(1 + pl^uo{9d)), 9s is the static contact angle. 



Uo[Ud) = 



sm9d-edCosed ^. /7pg(l-H4a^)\ 



1/2 



smUdCOsUd-fd V tI^ J 

is the characteristic speed associated with the parameters that the 'additional' physics of wetting brings in to 

resolve the moving contact-line problem, p^, pL, a, 0, 7, r are material constants characteriz ing the contacting 



2002 : 



Shikhmurzaev 



20071 ). The comparison 



media whose values can be found elsewhere (jBlake fc Shikhmurzaev 

of (12.4p with experimental data published in the literature has shown a very good agreement ([Shikhmurzaev 



20071 ) so that (|2.4p can be regarded as a reliable representation of the dynamic contact angle behaviour. 



Wetting front in a porous medium 5 

An implicit assumption we made above, namely that in the wetting mode the contact-line speed ui equals 

to the cross-sectionally averaged flow velocity uijiow associated with the motion of the meniscus as a whole, is 

not immediately obvious, especially given that, as the meniscus breaks through the threshold and goes into the 

wetting mode, initially, one has both the moving contact line and the varying shape of the meniscus. In other 

words, in general one has 

a dOd 

UlJlow = Ul + 7— . „ S2 -7-, 

(1 + smOdj at 
where, as before, we used that the meniscus has the shape of a spherical cap, albeit with a time-dependent 

radius of curvature. Formally, one can add this equation to the model together with an extra variable uijiow, 

which should replace ui in the text below, but this generalization would be beyond the accuracy of the model. 

Indeed, if L, U and T — L/U are the length, velocity and time scales for the macroscopic (Darcy-level) motion 

for which we are deriving the boundary conditions, then one has that the difference between uijiow and ui is 

of 0{a/L) and hence negligible in the continuum limit a/L — > 0. The possible deviation of the meniscus shape 

from a spherical cap takes place also on a vanishing scale. Therefore, within the accuracy of 0(1) as a/L -^ 0, 

in what follows we use that in the wetting mode uijiow = ui. It is worth pointing out here that the pores (i.e. 

capillaries) where the meniscus propagates in the wetting mode are assumed to be long compared to a, so that 

there is room for the meniscus to propagate in the wetting mode as it is described above. 

The speed ui at which individual menisci propagate in the wetting mode is not equal to the normal component 

of the velocity of the wetting front as a whole Un — n ■ u|goi (n is an outward normal) since the menisci 

intermittently go through both modes of motion and hence, on the Darcy scale, u„ must have contributions 

from both ui and the flow speed U2 associated with the threshold mode. Then, for m„ one has an equation 

u„ = Aiui + A2U2, (2.5) 

which is similar to (|2.ip . with the contribution of the rth mode proportional to its 'weight' Ai. In an isotropic 
porous medium, the 'weight' Ai of each mode of motion is essentially the relative time the meniscus spends 
in this mode. If Si is the fraction of the length on the pore scale corresponding to the ith mode of motion 
(si + S2 — 1), then the normalized time that the meniscus spends in this mode is Si/ui and hence, given that 
1/un = si/ui + S2/U2 and Ai + A2 = 1, one has 

A, = ^i"^ , A2 = ^^"^ ■ (2.6) 

S2U1 + S1U2 S2U1 + S1U2 

Then, as one would expect, the slowest (controlling) mode of motion tends to make a greater contribution to 

the pressure at the wetting front and the front's velocity. 
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Now, consider the threshold mode. When the moving meniscus runs into a barrier associated with the threshold 

mode, such as the edge at the end of a capillary or an asperity, the contact line gets pinned and the meniscus 

begins to deform until the contact angle reaches a certain value 0^, at which the contact line can move forward 

again (Fig. [ij. This is the essence of the threshold mode, and it comes into play only when 9d at which the 

meniscus arrives at a barrier is less than 0*. In other words, 

1, Od -e,>o 
si{ed,e,)^l , S2-1-S1, (2.7) 

[^ sio, 0d~O^ <0 
where sio (< 1) is a characteristic of the porous matrix. To find the functional dependence of parameters in the 
threshold mode, consider the dynamics of an individual meniscus. Assume that at a distance I upstream from 
the meniscus that has just run into a barrier and got its contact line pinned there is a (constant throughout the 
process) pressure p greater than —2a cos 9d/ a. Then, the meniscus, with the contact line which is now unable to 
move, will give in to this pressure and deform (Fig. [1]). We need to find the flow velocity U2 as an average over 
the cross-section and over the time required for the the contact angle to vary from 9d, at which the meniscus 
arrived at the barrier, to 0*, at which the contact line can advance again, and the pressure p2 as an average 
over the time of this process. 

Neglecting the contribution to the pressure drop due to the deviation of the velocity profile on the pore scale 
from parabolic in the immediate neighboorhood of the meniscus and calculating the flow rate in the capillary 
on an assumption that the meniscus retains the shape of a spherical cap (with a varying radius of curvature) 
throughout the process of its deformation, from the Stokes equation on the pore scale one has that the contact 

angle 9 satisfies an equation 

1 d9 a (pa 



a_ n 



cost 



(l + sin6')2 dt V V2ct 
which, if multiplied by a, essentially equates the fiow velocity averaged over a cross-section to the pressure 

gradient with a coefficient of proportionality corresponding to the parabolic profile in the pipe flow. Given that 

there is no discontinuity in the average flow velocity when the contact line becomes instantly pinned and the 

meniscus starts going from the wetting into the threshold mode, we have an equation 

aa /pa ^ \ 

_(^_+COS0,j=U„ 

which can be used to eliminate I. Now, for p2, i.e. the pressure —2(7COs9{t)/a averaged over a time interval 

T = aur^(g+cos0,)/(0,,0,;g 



(^'^'^-S 



d9 



g^ {I + sin 9)^ {pa / {2a) + COS 9)'- 
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which is needed for 9 to vary from 9d to 0*, one has 






where 

and [/]^ = f{b) — /(a). A similar procedure yields the velocity in the threshold mode averaged over T as 



1 /9 tt\ 1 ^/9 n 
- tan + - tan"^ 

2 V2 4/ 6 V2 4 



^^+,os9,)l[9,,9.;-) 
Now, the quantity p in equations (J2.8I) and (j2.9p must be specified in macroscopic terms. On the pore scale, 
it is the excess of p over the threshold capillary pressure p* = — 2a cos 9^,/ a that allows the meniscus to break 
through the barrier and go into the wetting regime again, with the contact line moving. If p reduces to p^ , then 
the meniscus comes to a halt, and, given that all menisci are modelled as the same, each of them will meet a 
similar barrier within a time negligible on the macroscopic scale and the wetting front as a whole will come to 
a stop. Formally, we have that, if p \ p*, then / — >■ oo and hence, according to (|2.8p and (|2.9I) . P2 — 5" P and 
M2 -^ 0. Then, we have from (|2.6p that Ai —^ 0, A2 -^ I, so that (|2.5p and (|2.ip yield that on the Darcy scale 
M„ — )■ and p ^^ p. Thus, macroscopically (i.e. on the Darcy scale) p is the stagnation pressure^ i.e. the pressure 
that one would have if the wetting front were at rest in its current position. In other words, macroscopically, 
we need to solve a conjugate problem 

V2p = 0, (rer!); n • Vp|aOi = 0, (2.10) 

with the boundary condition for p on dU,2 being the same as for p. Then, the value p|ani is the one we need to 
use in (|2.8p and (|2.9p . i.e. it is the pressure that builds on a meniscus whose contact line has been pinned. The 
conjugate problem must be solved in parallel with the main one as the latter requires the value of pj^Oi at the 
wetting front throughout its movement. 

Now, we can summarize the model as follows. In order to describe the propagation of the wetting front, one 
has to consider the bulk equations (II. ip and (|1.2p in the domain Vl with some boundary condition on 8^2 that 
specifies a particular problem, with the kinematic boundary condition at the wetting front 9f2i given by p.3p 
and the dynamic one by (j2.ip . The pressures pi and p2 that feature in (|2.ip are determined from (j2.2p . (|2.8p . 
with the coefficients Ai, A2 specified by (|2.6p . (12.71) . For the three variables 9d, ui and U2 appearing in ()2.2p . 
(|2.6p . (|2.7p and (|2.8I) one has three equations: (|2.3p (in particular (|2.4p ). (|2.5p and (|2.9p . whereas the pressure 
p\dQ,i featuring in (|2.8p and (|2.9p has to be found from the conjugate problem (|2.10p with the same boundary 
condition for p on 8^2 as for p. If gravity is to be taken into account, one has to replace p and p in (jl.ip and 



8 Y.D. Shikhmurzaev and J.E. Sprittles 

(|2.10p with {p + pgz) and {p + pgz), respectively (p is the fluid's density, g is the gravitational acceleration, 

z is the coordinate directed against gravity). Besides the bulk permeability k, the geometry of the porous 

matrix enters the model via three effective parameters, sio, a and 6^. The above model can be generalized by 

incorporating threshold modes associated with different values of 9^, , replacing the step- function (|2.7p for each 

of them by a more complex one to account for the types of thresholds and the corresponding generalization of 

(j2.6p to incorporate various possibilities for the pore network topology. 

3. An illustrative example 

Consider the unsteady one-dimensional imbibition against gravity, with z = h[t) as the position of the wetting 

front [dU,i ) , p and p in (|l.ll) and (I2.10p replaced with p + pgz and p + pgz respectively, and p = p = pa dX z = Q 

as a boundary condition on d^2- Then, Laplace's equations (jl.2p and (I2.10p for p and p in the one-dimensional 

case yield that both p and p are linear functions of z, and from the conjugate problem (j2.10p with the condition 

p — Pq a,i z = Q wc have that p{z, t) ~ p^ ~ pgz. Then, using for p its linear dependence on z and the same 

condition on d^2i i-e. p = p^ ai z = 0, we can express the pressure gradient dp/dz at z = /i in terms of the 

current position of the wetting front: dp{h,t)/dz = {p{h,t) — p^)/h. Using this expression in Darcy's equation 

(jl.ip . where, as mentioned above, p is replaced with p + pgz, and substituting the latter into the kinematic 

condition p.3p . which now can be written down simply as dh/dt — Un{h,t), we arrive at 

dh K f pq — p(h,t) \ 

dt p.\ h J 

This equation together with algebraic equations (12.11) . where p\dni = p{h,t), (|2.2p . (|2.4p . (|2.5p . where u„ = 

dh/dt, (|2.6I) . (|2.7p . (|2.8p . where we now use the solution of the conjugate problem p\dni = Po — pgh, and 

(|2.9p form a closed system for h, p{h,t), pi, 9d, ui, Ai, A2, si, S2, P2 and U2. Typical curves representing the 

dependence of h, scaled with Hq — 2a/{pga) and rising from the initial position h{0) = 0, on time t, scaled with 

To = 2ap/((pg)'^aK), are shown in Fig. [2] with the v alues of paramete rs given in the figure caption. Comparing 



these dependencies (curves 2-5) with Washburn's (jWashburn 



19211 ) curves for a meniscus propagating with 



9d = 9s in a capillary (curve 0: no gravity; curve 1: gravity included), we can see that the present model, besides 
realistically giving a finite speed of rise at the onset of imbibition, predicts a slightly lower rate of propagation 
of the wetting front as it is slowed down by both the velocity-dependence of the contact angle and the presence 
of the threshold mode. The latter comes into play once 9d < 9^, and dictates that the wetting front will come 
to a halt before it reaches the maximum possible height of imbibition h^ax — 2a cos 9s/ (pga) and, importantly, 
contrary to how curves 1 and 2 approach h^ax asymptotically, this coming to a halt occurs in a finite time. 
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Figure 2. Time-dependence of ID imbibition calculated using the derived model. 0: Washburn, no gravity; 1: Washburn, 
gravity included; 2: sw = 1; 3: sio = 0.1, 9* = 30°; 4: sio = 0.9, 6, = 60°; 5: sio = 0.1, 61* = 60°. For all curves 6*^ = 0°, 
po = 0, and fiUd / (upg) = 10^, pl^ = 0.6 for curves 2-5. 

In this position, which is determined only by 6^ (see curves 4 and 5), the wetting front still has the capacity 
to propagate provided that some other physical mechanism helps it to break through the threshold. In this 
connection, it is worth pointing out that, as has been observed experimentally, for some systems, there is a 
change of regime from an essentially Washburn-t ype to a completely different one halfway between the onset of 



the process and the maximum imbibition height (jDelker et al. 



1996; 



Lago fc Arauio 



20011 ). so that the present 



example, considered here as an illustration of how the developed approach works, provides a 'building block' 
for the modelling of this, yet unexplained, phenomenon. 

4. Concluding remarks 

The developed approach and, in particular, the technique of conjugate problems it uses to incorporate the 
threshold mode of motion provide a transparent framework for building models of increased complexity. The 
threshold mode is the key to describing such effects as formation of trapped pockets of the displaced fluid 
that can be left behind the wetting front and their subsequent dynamics without resorting to thermodynamic 
arguments and the necessity to specify, increasingly multi-parametric, thermodynamic dependencies for this, 
basically mechanical, system. It should be noted, however, that the next step in the developing of the model 
towards greater complexity, sketched at the end of Section[2j is nontrivial as, for different values of the threshold 
angle 0*, it becomes necessary to bring in and specify the topology of the porous matrix with respect to the 
connectivity of the wetting front. 
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